% Estimate Bid-Ask illiquidity spread

% A -- is my bid_ask_matrix of bid-ask prices 
% odd columns == bid prices
% even columns == ask prices

% Get the number of columns in A
[rows, cols] = size(A);

% Extract all the odd-numbered columns == Bid prices
bid = A(:, 1:2:cols);

% Extract all the even-numbered columns == Ask prices
ask = A(:, 2:2:cols);

% Estimate bid ask-spread

bid_ask_spread= (ask - bid) ./ ask;  % Calculate the spread using (Ask - Bid) / Ask

% Reshape the bid_ask_spread into single column vector 

liquidity=bid_ask_spread(:); 
%_________________________________________________________________________


% Get the log returns of CDS 

r= diff(log(cds));
% Get the Squared returns

R=r.^2;

% Get the RV of squared log returns for each month (120) 
%_________________________________________________________________________
y = 1;
for i = 0:21.71666:2606
RV_m(y,: ) = sum(R(i+1:i+21.7166,:));  
y=y+1;
end

% Reshape RV into single column vector
RV=RV_m(:);

% Get the Bipower Variation BV sum (from 2 to M)
%__________________________________________________________________________
rr=abs(R); % Get the absolute log returns
p1=rr; 
y=1;
for i=0:1:99999999
    C1(y,:)=prod(p1(i+1:i+2,:)); % I need C1 only for to calculate C - drop C1 
    y=y+1;
end

p2=C1;
y=1;
%global C 
for i = 0:21.70833:999999999
C(y,:) = sum(p2(i+2:i+21.70833,:)); % this is the sum (from 2 to M)that I need for BV
y=y+1;
end

% BV1 = ComputeBV1;
%__________________________________________________________________________
M = 21.70833; % M = 9; % number of intraday (intramonth) observations (for one day, not the whole sample)
delta = 1./M;
RR = r.^(4/3); % the abs(return).^4/3 see also TP appendix A
EE=RR;


mu1= (pi/2); % eq 4 or see Appendix A

BV1 = mu1.^(-2) .* C;

% TP = Compute Tripower Quarticity TQ sum (from 3 to M)
%__________________________________________________________________________

p3=R;
y=1;
for i=0:1:18832
    E1(y,:)=prod(p3(i+1:i+3,:)); % I need E1 only for to calculate E- drop E1 
    y=y+1;
end

p4=E1;
y=1;

for i = 0:15.166:18832
EE(y,:) = sum(p4(i+3:i+15.166,:)); % this is the sum (from 3 to M)that I need for TQ (TP)
y=y+1;
end
mu = 1./(4.*delta.*((gamma(7/6))./(gamma(1/2))).^3); % see appendix A

TP = mu.*EE; % E or EE is the: R = r.^(4/3); the abs(return).^4/3 see also TP appendix A

% RJ = Compute Realized Jumps 
%_________________________________________________________________________
 RJ = (RV-BV1)./RV;

% Z1 = Compute the test statistics 

% Zhang, Zhou, Zhu (RFS, 2009), Appendix A 
% function Z1 = ComputeZ1(~, RV, BV1, ~, TP)
% Copmutes the test statistics for assessing the significance of the daily
% jump component. Thus, only statistically extreme positive values of RV-BV/RV
% are attributed to the jump component. BNS(2004), Huang and Tauchen(2005)
% and Andersen, Bollerslev and Diebold(2007) for significant jumps test
% statistics based in ratio statistics
%__________________________________________________________________________

Z1 = RJ./(((pi/2).^2 + pi - 5)* delta * max(1, TP./(BV1.^2) )).^(1/2); % Z1 

% RV is the Realized Volatility
% BV1 is the Bipower Variation
% TP is the Tripower Quarticity (practical proxy for the Integrated Variance)

% Realized Jump Volatility == Computes the stat significance of Jumps
%__________________________________________________________________________
% 95% significance level
for i = 1:1:size(Z1,1)
    for j = 1:1:size(Z1,2)
        if Z1(i,j) >  1.644854 %1.644854  is of 95% 3.090232 %99.9%[the correct] % 2.326348 is 99% significance level
RVJ(i,j) = sqrt(RV(i,j))-sqrt(BV1(i,j));
RVC(i,j) = sqrt(BV1(i,j));
        end
    if Z1(i,j) < 1.644854 
   RVJ(i,j) = 0;
    RVC(i,j) = sqrt(RV(i,j));     
    end
    end
end

% 99 % Significance level
for i = 1:1:size(Z1,1)
    for j = 1:1:size(Z1,2)
        if Z1(i,j) >  2.326348 % 2.326348 is 99% significance level
RVJ_99(i,j) = sqrt(RV(i,j))-sqrt(BV1(i,j));
RVC_99(i,j) = sqrt(BV1(i,j));
        end
    if Z1(i,j) < 2.326348 
   RVJ_99(i,j) = 0;
    RVC_99(i,j) = sqrt(RV(i,j));     
    end
    end
end